Preparation_Amazon_data
2022-12-28
1 Load libraries
# load libraries
library(tictoc)
library(tidyverse)
library(terra)
library(sf)
library(RColorBrewer)
library(rts)
library(doParallel)
library(parallel)
library(foreach)
library(data.table)2 Functions
# Function to save layers
saveLayers = function(dataRast, shapeFile, pathFile, strRemove,
x_min=xmin(dataRast), x_max=xmax(dataRast),
y_min=ymin(dataRast), y_max=ymax(dataRast)){
# Iterate all the layers
for (i in 1:nlyr(dataRast)){
print(i)
dataRast.i <- dataRast[[i]]
# Crop data
dataRast.basin.i <- mask(dataRast.i, mask = shapeFile) %>%
crop(ext(x_min, x_max, y_min, y_max))
# Path of the file
path <- gsub(strRemove, "", names(dataRast.basin.i)) %>%
paste0(pathFile, ., ".tif")
# Save file
writeRaster(dataRast.basin.i, path, overwrite=TRUE)
# Clear memory
gc()
}
}
# Function to rename layers
renameLayers <- function(dataRast, fileStart, prefix){
# index of file name in the path
id.date <- unlist(gregexpr(fileStart, sources(dataRast)[[1]])) + nchar(fileStart)
# rename layers
names(dataRast) <-
substr(sources(dataRast), id.date, (id.date+50)) %>%
gsub(".tif", "_1", .) %>%
as.Date("%Y_%m_%d") %>%
format(., '%Y_%m') %>%
paste0(prefix, .)
return(dataRast)
}3 Load data
load("R_data/data-amaz_na3.RData")#=====================
4 Import shape file
# Import shape file
amaz.basin.shp <- st_read("Amazon_new_data/0. Amazon_shapefile/projected/amazon_shp_projected.shp")## Reading layer `amazon_shp_projected' from data source
## `/Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/0. Amazon_shapefile/projected/amazon_shp_projected.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -2156811 ymin: 1625314 xmax: 1745999 ymax: 4555427
## Projected CRS: South_America_Albers_Equal_Area_Conic#=====================
5 Burnt Area data
5.1 Import data
# list of files
amaz.burntArea.list <- list.files("Amazon_new_data/1. Burnt Area/03. Working Data",
full.names=TRUE,
pattern = ".tif$")
# Import data with "Terra"
burntArea.rast <- rast(amaz.burntArea.list)
burntArea.rast## class : SpatRaster
## dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## sources : burntarea_working_2001_1.tif
## burntarea_working_2001_10.tif
## burntarea_working_2001_11.tif
## ... and 237 more source(s)
## names : fire_~_proj, fire_~_proj, fire_~_proj, fire_~_proj, fire_~_proj, fire_~_proj, ...
## min values : -2, -2, -2, -2, -2, -2, ...
## max values : 1, 1, 1, 1, 1, 1, ...5.2 Verification of the data
5.2.1 Names of layers
names(burntArea.rast) %>% head()## [1] "fire_2001_05_06_1_proj" "fire_2001_05_06_10_proj"
## [3] "fire_2001_05_06_11_proj" "fire_2001_05_06_12_proj"
## [5] "fire_2001_05_06_2_proj" "fire_2001_05_06_3_proj"5.2.2 Find duplicate name layers
idx.dplc <- names(burntArea.rast) %>% duplicated() %>% which()
idx.dplc## [1] 1445.2.3 Print duplicate name layers
burntArea.rast[[names(burntArea.rast) == names(burntArea.rast)[idx.dplc]]]## class : SpatRaster
## dimensions : 5860, 7806, 2 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## sources : burntarea_working_2012_7.tif
## burntarea_working_2012_9.tif
## names : fire_2012_05_06_9_proj, fire_2012_05_06_9_proj
## min values : -1, -1
## max values : -1, -15.2.4 Verification of the values
# Verification of the values
burntArea.minmax <- minmax(burntArea.rast) %>% t() %>% as.data.frame()
burntArea.minmaxburntArea.minmax[which((burntArea.minmax[,1] != -2) & (burntArea.minmax[,2] != 1)),]5.2.5 Plot of the months September and November 2012
# Create palette of color
pal <- colorRampPalette(c("lightblue", "forestgreen", "firebrick"))
# burntArea.freq <- freq(burntArea.rast, digits=0, usenames=T)
burntArea.freq[(burntArea.freq[,2] != -2)&(burntArea.freq[,2] != 0)&(burntArea.freq[,2] != 1),]plot(burntArea.rast[[c(142,144)]], col=pal(3))
5.3 Modification of the data: Rename layers
5.3.1 Rename layers
# Rename layers
burntArea.rast <- renameLayers(burntArea.rast, 'burntarea_working_', 'fire_')
# Order layers
ordered.names.1 <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month") %>%
format(., '%Y_%m') %>%
paste0("fire_", .)
burntArea.rast <- burntArea.rast[[ordered.names.1]]
burntArea.rast## class : SpatRaster
## dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## sources : burntarea_working_2001_1.tif
## burntarea_working_2001_2.tif
## burntarea_working_2001_3.tif
## ... and 237 more source(s)
## names : fire_2001_01, fire_2001_02, fire_2001_03, fire_2001_04, fire_2001_05, fire_2001_06, ...
## min values : -2, -2, -2, -2, -2, -2, ...
## max values : 1, 1, 1, 1, 1, 1, ...5.3.2 Create ‘rts’ object
# Create a sequence date for 'rts' object
seq.dates <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month")
burntArea.rts <- rts(burntArea.rast, seq.dates)
burntArea.rts## Raster Time Series with monthly periodicity from 2001-01-01 to 2020-12-01
## class : SpatRasterTS
## raster filename : /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/1. Burnt Area/03. Working Data/burntarea_working_2001_1.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/1. Burnt Area/03. Working Data/burntarea_working_2001_2.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/1. Burnt Area/03. Working Data/burntarea_working_2001_3.tif ...
## raster dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## raster resolution : 500, 500 (x, y)
## raster extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## min values : -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 ...
## max values : 1 1 1 1 1 1 1 1 1 1 ...5.4 Plot
5.4.1 Define palette of color
# Create palette of color
pal <- colorRampPalette(c("lightblue", "forestgreen", "firebrick"))
pal2 <- colorRampPalette(c("forestgreen", "firebrick"))5.4.2 Plot all the area
# Plot Amazon in the year 2012.
plot(burntArea.rts[['2012']], col=pal(3))
# plot Amazon in October 2020.
plot(burntArea.rts[['2020-10-01']], col=pal(3), main="Burnt area in October 2020")
plot(amaz.basin.shp$geometry, add=T)
#
burntArea.rts[['2020-10-01']] %>%
mask(mask = amaz.basin.shp) %>%
plot(col=pal(3), main="Burnt area in October 2020")##
|---------|---------|---------|---------|
=========================================
plot(amaz.basin.shp$geometry, add=T)
5.4.3 Plot south of Amazon
x_min=-1.25e+06; x_max=0.05e+06; y_min=ymin(burntArea.rast); y_max=2.4e+06
# Plot south of Amazon in the year 2020.
burntArea.rts[['2020']]@raster %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
mask(mask = amaz.basin.shp) %>%
plot(col=pal(3))##
|---------|---------|---------|---------|
=========================================
|---------|---------|---------|---------|
=========================================
# Plot south of Amazon in October 2020.
burntArea.rts[['2020-10-01']] %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
mask(mask = amaz.basin.shp) %>%
plot(col=pal(3), main="Burnt area in October 2020")
plot(amaz.basin.shp$geometry, add=T)
5.5 Repalce undifined values by NA
burntArea.rast[[c(139, 141)]][burntArea.rast[[c(139, 141)]] == -1] <- NA##
|---------|---------|---------|---------|
=========================================
|---------|---------|---------|---------|
=========================================
burntArea.rast[[c(139, 141)]]## class : SpatRaster
## dimensions : 5860, 7806, 2 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## sources : spat_r1U8IiItX34kjKM_10135.tif
## spat_r1U8IiItX34kjKM_10135.tif
## names : fire_2012_07, fire_2012_09
## min values : NaN, NaN
## max values : NaN, NaN5.6 Verification of the NA values
# cl <- makeCluster(detectCores() - 1)
# registerDoParallel(cl, cores=detectCores() - 1)
#
# tic("burntArea")
# #Raster to datatable in parallel: one raster per thread
# rasList <- foreach (ras_id=amaz.burntArea.list, .packages=c('terra', 'sf'), .combine='c') %dopar% {
# #Read all rasters into one big stack
# ras <- rast(ras_id)
# # Rename layers
# ras <- renameLayers(ras, 'burntarea_working_', '')
# # Replace negative value by `NA`
# if (names(ras) %in% c("2012_07", "2012_09")) {ras[ras == -1] <- NA}
# #
# ras.nonNA <- not.na(ras)
# ras.nonNA.mask <- mask(ras.nonNA, amaz.basin.shp)
# ras.freq.na <- freq(ras.nonNA.mask, digits=0, value=0, usenames=T)
#
# list(ras.freq.na)
# }
# stopCluster(cl)
# toc()
#
# #Bind all per-raster into one dataframe
# burntArea.freq.na <- rbindlist(rasList, fill=T, use.names=T)
# #
# colnames(burntArea.freq.na)[3] <- "burntArea_na"
# burntArea.freq.na <- burntArea.freq.na[order(burntArea.freq.na$layer)]
burntArea.freq.na#=====================
6 Land Cover data
6.1 Import data
# list of files
amaz.landCover.list <- list.files("Amazon_new_data/2. Land Cover/03. Working Data",
full.names=TRUE,
pattern = ".tif$")
# Import data with "Terra"
landCover.rast <- rast(amaz.landCover.list)
landCover.rast## class : SpatRaster
## dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## sources : landcover_working_2001_1.tif
## landcover_working_2001_10.tif
## landcover_working_2001_11.tif
## ... and 237 more source(s)
## names : lulc_~_proj, lulc_~_proj, lulc_~_proj, lulc_~_proj, lulc_~_proj, lulc_~_proj, ...
## min values : 0, 0, 0, 0, 0, 0, ...
## max values : 10, 10, 10, 10, 10, 10, ...6.2 Verification of the data
6.2.1 Verification the Names of layers
names(landCover.rast) %>% head(24)## [1] "lulc_2001_proj" "lulc_2001_proj" "lulc_2001_proj" "lulc_2001_proj"
## [5] "lulc_2001_proj" "lulc_2001_proj" "lulc_2001_proj" "lulc_2001_proj"
## [9] "lulc_2001_proj" "lulc_2001_proj" "lulc_2001_proj" "lulc_2001_proj"
## [13] "lulc_2002_proj" "lulc_2002_proj" "lulc_2002_proj" "lulc_2002_proj"
## [17] "lulc_2002_proj" "lulc_2002_proj" "lulc_2002_proj" "lulc_2002_proj"
## [21] "lulc_2002_proj" "lulc_2002_proj" "lulc_2002_proj" "lulc_2002_proj"6.2.2 Verification of the values
# Verification of the values
landCover.minmax <- minmax(landCover.rast) %>% t() %>% as.data.frame()
landCover.minmax6.2.3 Frequency of each value
# landCover.freq <- freq(landCover.rast, digits=0, usenames=T)
landCover.freq#
landCover.freq[which(!(landCover.freq$value %in% c(0,1,2,3,4,5,6,7,8,9,10)) ),]6.3 Modification of the data
6.3.1 Rename layers
# Rename layers
landCover.rast <- renameLayers(landCover.rast, 'landcover_working_', 'lulc_')
# Order layers
ordered.names.2 <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month") %>%
format(., '%Y_%m') %>%
paste0("lulc_", .)
landCover.rast <- landCover.rast[[ordered.names.2]]
landCover.rast## class : SpatRaster
## dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## sources : landcover_working_2001_1.tif
## landcover_working_2001_2.tif
## landcover_working_2001_3.tif
## ... and 237 more source(s)
## names : lulc_2001_01, lulc_2001_02, lulc_2001_03, lulc_2001_04, lulc_2001_05, lulc_2001_06, ...
## min values : 0, 0, 0, 0, 0, 0, ...
## max values : 10, 10, 10, 10, 10, 10, ...6.3.2 Create ‘rts’ object
# Create a sequence date for 'rts' object
seq.dates <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month")
landCover.rts <- rts(landCover.rast, seq.dates)
landCover.rts## Raster Time Series with monthly periodicity from 2001-01-01 to 2020-12-01
## class : SpatRasterTS
## raster filename : /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/2. Land Cover/03. Working Data/landcover_working_2001_1.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/2. Land Cover/03. Working Data/landcover_working_2001_2.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/2. Land Cover/03. Working Data/landcover_working_2001_3.tif ...
## raster dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## raster resolution : 500, 500 (x, y)
## raster extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## min values : 0 0 0 0 0 0 0 0 0 0 ...
## max values : 10 10 10 10 10 10 10 10 10 10 ...6.4 Plot
6.4.1 Plot all the area
# Plot Amazon in the year 2020
plot(landCover.rts[['2020']], col=colorRamps::blue2green(10))
# plot Amazon in October 2020.
plot(landCover.rts[['2020-10-01']], col=colorRamps::blue2green(10), main="Land cover in October 2020")
plot(amaz.basin.shp$geometry, add=T)
#
landCover.rts[['2020-10-01']] %>%
mask(mask = amaz.basin.shp) %>%
plot(col=colorRamps::blue2green(10), main="Land cover in October 2020")##
|---------|---------|---------|---------|
=========================================
plot(amaz.basin.shp$geometry, add=T)
6.4.2 Plot south of Amazon
x_min=-1.25e+06; x_max=0.05e+06; y_min=ymin(burntArea.rast); y_max=2.4e+06
# Plot south of Amazon in the year 2020.
landCover.rts[['2020']]@raster %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
mask(mask = amaz.basin.shp) %>%
plot(col=colorRamps::blue2green(10))##
|---------|---------|---------|---------|
=========================================
|---------|---------|---------|---------|
=========================================
# Plot south of Amazon in October 2020.
landCover.rts[['2020-10-01']] %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
mask(mask = amaz.basin.shp) %>%
plot(col=colorRamps::blue2green(10), main="Land cover in October 2020")
plot(amaz.basin.shp$geometry, add=T)
6.5 Verification of the NA values
# cl <- makeCluster(detectCores() - 1)
# registerDoParallel(cl, cores=detectCores() - 1)
#
# tic("landCover")
# #Raster to datatable in parallel: one raster per thread
# rasList <- foreach (ras_id=amaz.landCover.list, .packages=c('terra', 'sf'), .combine='c') %dopar% {
# #Read all rasters into one big stack
# ras <- rast(ras_id)
# # Rename layers
# ras <- renameLayers(ras, 'landcover_working_', '')
# #
# ras.nonNA <- not.na(ras)
# ras.nonNA.mask <- mask(ras.nonNA, amaz.basin.shp)
# ras.freq.na <- freq(ras.nonNA.mask, digits=0, value=0, usenames=T)
#
# list(ras.freq.na)
# }
# stopCluster(cl)
# toc()
#
# #Bind all per-raster into one dataframe
# landCover.freq.na <- rbindlist(rasList, fill=T, use.names=T)
# #
# colnames(landCover.freq.na)[3] <- "landCover_na"
# landCover.freq.na <- landCover.freq.na[order(landCover.freq.na$layer)]
landCover.freq.na#=====================
7 Precipitation data
7.1 Import data
# list of files
amaz.precipitation.list <- list.files("Amazon_new_data/3. Precipitation/03. Working Data",
full.names=TRUE,
pattern = ".tif$")
# Import data with "Terra"
precipitation.rast <- rast(amaz.precipitation.list)
precipitation.rast## class : SpatRaster
## dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## sources : precipitation_working_2001_1.tif
## precipitation_working_2001_10.tif
## precipitation_working_2001_11.tif
## ... and 237 more source(s)
## names : prec_~_proj, prec_~_proj, prec_~_proj, prec_~_proj, prec_~_proj, prec_~_proj, ...
## min values : 0, 0, 0, 0, 0, 0, ...
## max values : 1461, 1735, 1828, 1934, 1433, 1390, ...7.2 Verification of the data
7.2.1 Verification the Names of layers
names(precipitation.rast) %>% head()## [1] "prec_2001_1_proj" "prec_2001_10_proj" "prec_2001_11_proj"
## [4] "prec_2001_12_proj" "prec_2001_2_proj" "prec_2001_3_proj"7.2.2 Verification of the values
# Verification of the values
precipitation.minmax <- minmax(precipitation.rast) %>% t() %>% as.data.frame()
precipitation.minmax7.3 Modification of the data
7.3.1 Rename layers
# Rename layers
precipitation.rast <- renameLayers(precipitation.rast, 'precipitation_working_', 'prec_')
# Order layers
ordered.names.3 <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month") %>%
format(., '%Y_%m') %>%
paste0("prec_", .)
precipitation.rast <- precipitation.rast[[ordered.names.3]]
precipitation.rast## class : SpatRaster
## dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## sources : precipitation_working_2001_1.tif
## precipitation_working_2001_2.tif
## precipitation_working_2001_3.tif
## ... and 237 more source(s)
## names : prec_2001_01, prec_2001_02, prec_2001_03, prec_2001_04, prec_2001_05, prec_2001_06, ...
## min values : 0, 0, 0, 0, 0, 0, ...
## max values : 1461, 1433, 1390, 1096, 1678, 1706, ...7.3.2 Create ‘rts’ object
# Create a sequence date for 'rts' object
seq.dates <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month")
precipitation.rts <- rts(precipitation.rast, seq.dates)
precipitation.rts## Raster Time Series with monthly periodicity from 2001-01-01 to 2020-12-01
## class : SpatRasterTS
## raster filename : /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/3. Precipitation/03. Working Data/precipitation_working_2001_1.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/3. Precipitation/03. Working Data/precipitation_working_2001_2.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/3. Precipitation/03. Working Data/precipitation_working_2001_3.tif ...
## raster dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## raster resolution : 500, 500 (x, y)
## raster extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## min values : 0 0 0 0 0 0 0 0 0 0 ...
## max values : 1461 1433 1390 1096 1678 1706 1807 1988 2230 1735 ...7.4 Plot
7.4.1 Plot all the area
# Plot Amazon in the year 2020
plot(precipitation.rts[['2020']], col=grDevices::cm.colors(1000))
# plot Amazon in October 2020.
plot(precipitation.rts[['2020-10-01']], col=grDevices::cm.colors(1000), main="Precipitation in October 2020")
plot(amaz.basin.shp$geometry, add=T)
#
precipitation.rts[['2020-10-01']] %>%
mask(mask = amaz.basin.shp) %>%
plot(col=grDevices::cm.colors(1000), main="Precipitation in October 2020")##
|---------|---------|---------|---------|
=========================================
plot(amaz.basin.shp$geometry, add=T)
7.4.2 Plot south of Amazon
x_min=-1.25e+06; x_max=0.05e+06; y_min=ymin(burntArea.rast); y_max=2.4e+06
# Plot south of Amazon in the year 2020.
precipitation.rts[['2020']]@raster %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
mask(mask = amaz.basin.shp) %>%
plot(col=grDevices::cm.colors(1000))##
|---------|---------|---------|---------|
=========================================
|---------|---------|---------|---------|
=========================================
# Plot south of Amazon in October 2020.
precipitation.rts[['2020-10-01']] %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
mask(mask = amaz.basin.shp) %>%
plot(col=grDevices::cm.colors(1000), main="Precipitation in October 2020")
plot(amaz.basin.shp$geometry, add=T)
7.5 Verification of the NA values
# cl <- makeCluster(detectCores() - 1)
# registerDoParallel(cl, cores=detectCores() - 1)
#
# tic("precipitation")
# #Raster to datatable in parallel: one raster per thread
# rasList <- foreach (ras_id=amaz.precipitation.list, .packages=c('terra', 'sf'), .combine='c') %dopar% {
# #Read all rasters into one big stack
# ras <- rast(ras_id)
# # Rename layers
# ras <- renameLayers(ras, 'precipitation_working_', '')
# #
# ras.nonNA <- not.na(ras)
# ras.nonNA.mask <- mask(ras.nonNA, amaz.basin.shp)
# ras.freq.na <- freq(ras.nonNA.mask, digits=0, value=0, usenames=T)
#
# list(ras.freq.na)
# }
# stopCluster(cl)
# toc()
#
# #Bind all per-raster into one dataframe
# precipitation.freq.na <- rbindlist(rasList, fill=T, use.names=T)
# #
# colnames(precipitation.freq.na)[3] <- "precipitation_na"
# precipitation.freq.na <- precipitation.freq.na[order(precipitation.freq.na$layer)]
precipitation.freq.na7.6 Sample of NA values (February 2002)
7.6.1 Compute
precipitation.2002.02.name <- c("Amazon_new_data/3. Precipitation/03. Working Data/precipitation_working_2002_2.tif")
precipitation.2002.02.rast <- rast(precipitation.2002.02.name)
precipitation.2002.02.rast## class : SpatRaster
## dimensions : 5860, 7806, 1 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## source : precipitation_working_2002_2.tif
## name : prec_2002_2_proj
## min value : 0
## max value : 1439#
precipitation.2002.02.nonNA <- not.na(precipitation.2002.02.rast)##
|---------|---------|---------|---------|
=========================================
precipitation.2002.02.nonNA.mask <- mask(precipitation.2002.02.nonNA, amaz.basin.shp)##
|---------|---------|---------|---------|
=========================================
#
precipitation.2002.02.nonNA.df <- as.data.frame(precipitation.2002.02.nonNA.mask, xy=T)
precipitation.2002.02.NA.df <- precipitation.2002.02.nonNA.df[precipitation.2002.02.nonNA.df[,3]==F,]
precipitation.2002.02.NA.dfdim(precipitation.2002.02.NA.df)[1]## [1] 347.6.2 Plot
# Plot
plot(precipitation.2002.02.rast)
plot(precipitation.2002.02.nonNA.mask, col=c("darkorange1", "forestgreen"))
plot(amaz.basin.shp$geometry, add=T)
#
x_min=-0.45e+06; x_max=-0.3e+06; y_min=ymin(precipitation.2002.02.nonNA.mask); y_max=1.71e+06
precipitation.2002.02.nonNA.mask %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
plot(col=c("darkorange1", "forestgreen"))
#=====================
8 Soil Moisture data
8.1 Import data
# list of files
amaz.soilMoisture.list <- list.files("Amazon_new_data/4. Soil Moisture/03. Working Data",
full.names=TRUE,
pattern = ".tif$")
# Import data with "Terra"
soilMoisture.rast <- rast(amaz.soilMoisture.list)
soilMoisture.rast## class : SpatRaster
## dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## sources : soilmoisture_working_2001_1.tif
## soilmoisture_working_2001_10.tif
## soilmoisture_working_2001_11.tif
## ... and 237 more source(s)
## names : soilm~_proj, soilm~_proj, soilm~_proj, soilm~_proj, soilm~_proj, soilm~_proj, ...
## min values : -9.990000e+08, -9.990000e+08, -9.990000e+08, -9.990000e+08, -9.990000e+08, -9.990000e+08, ...
## max values : 1.076058e+03, 4.290171e+03, 7.905507e+02, 5.636223e+02, 7.509722e+02, 6.363759e+02, ...8.2 Verification of the data
8.2.1 Verification the Names of layers
names(soilMoisture.rast) %>% head()## [1] "soilm_2001_1_proj" "soilm_2001_10_proj" "soilm_2001_11_proj"
## [4] "soilm_2001_12_proj" "soilm_2001_2_proj" "soilm_2001_3_proj"8.2.2 Verification of the values
# Verification of the values
soilMoisture.minmax <- minmax(soilMoisture.rast) %>% t() %>% as.data.frame()
soilMoisture.minmax8.2.3 Frequency
# soilMoisture.freq <- freq(soilMoisture.rast, digits=3, usenames=T)
soilMoisture.freq#
soilMoisture.freq[soilMoisture.freq$value < 0,]8.3 Modification of the data
8.3.1 Rename layers
# Rename layers
soilMoisture.rast <- renameLayers(soilMoisture.rast, 'soilmoisture_working_', 'soilm_')
# Order layers
ordered.names.4 <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month") %>%
format(., '%Y_%m') %>%
paste0("soilm_", .)
soilMoisture.rast <- soilMoisture.rast[[ordered.names.4]]
# Print
soilMoisture.rast## class : SpatRaster
## dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## sources : soilmoisture_working_2001_1.tif
## soilmoisture_working_2001_2.tif
## soilmoisture_working_2001_3.tif
## ... and 237 more source(s)
## names : soilm~01_01, soilm~01_02, soilm~01_03, soilm~01_04, soilm~01_05, soilm~01_06, ...
## min values : -9.990000e+08, -9.990000e+08, -9.990000e+08, -9.990000e+08, -9.990000e+08, -9.990000e+08, ...
## max values : 1.076058e+03, 7.509722e+02, 6.363759e+02, 7.215025e+02, 6.892701e+02, 5.409775e+02, ...cat('\n---------- \n\n')##
## ----------names(soilMoisture.rast) %>% head()## [1] "soilm_2001_01" "soilm_2001_02" "soilm_2001_03" "soilm_2001_04"
## [5] "soilm_2001_05" "soilm_2001_06"8.3.2 Create ‘rts’ object
# Create a sequence date for 'rts' object
seq.dates <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month")
soilMoisture.rts <- rts(soilMoisture.rast, seq.dates)
soilMoisture.rts## Raster Time Series with monthly periodicity from 2001-01-01 to 2020-12-01
## class : SpatRasterTS
## raster filename : /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/4. Soil Moisture/03. Working Data/soilmoisture_working_2001_1.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/4. Soil Moisture/03. Working Data/soilmoisture_working_2001_2.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/4. Soil Moisture/03. Working Data/soilmoisture_working_2001_3.tif ...
## raster dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## raster resolution : 500, 500 (x, y)
## raster extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## min values : -1e+09 -1e+09 -1e+09 -1e+09 -1e+09 -1e+09 -1e+09 -1e+09 -1e+09 -1e+09 ...
## max values : 1076 751 636 722 689 541 713 1167 1568 4290 ...8.4 Plot
8.4.1 Plot all the area
# Plot Amazon in the year 2020
plot(soilMoisture.rts[['2020']], col=colorRamps::matlab.like2(1000))
# plot Amazon in October 2020.
plot(soilMoisture.rts[['2020-10-01']], col=colorRamps::matlab.like2(1000), main="Soil Moisture in October 2020")
plot(amaz.basin.shp$geometry, add=T)
#
soilMoisture.rts[['2020-10-01']] %>%
mask(mask = amaz.basin.shp) %>%
plot(col=colorRamps::matlab.like2(1000), main="Soil Moisture in October 2020")##
|---------|---------|---------|---------|
=========================================
plot(amaz.basin.shp$geometry, add=T)
#
x_min=0.1e+06; x_max=0.5e+06; y_min=4.25e+06; y_max=ymax(burntArea.rast)
soilMoisture.rts[['2020-10-01']] %>%
mask(mask = amaz.basin.shp) %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
plot(col=colorRamps::matlab.like2(1000), main="Soil Moisture in October 2020")##
|---------|---------|---------|---------|
=========================================
plot(amaz.basin.shp$geometry, add=T)
8.4.2 Repalce negative values by NA for 2020
soilMoisture.2020.rts <- soilMoisture.rts[['2020']]
# Remove negative values.
soilMoisture.2020.rts@raster[soilMoisture.2020.rts@raster < 0] <- NA##
|---------|---------|---------|---------|
=========================================
|---------|---------|---------|---------|
=========================================
soilMoisture.2020.rts## Raster Time Series with monthly periodicity from 2020-01-01 to 2020-12-01
## class : SpatRasterTS
## raster filename : /private/var/folders/91/ckb_ktpx5mzcsls1l4tmm5xm0000gp/T/Rtmp6sflAg/spat_CKF46jJIMXAcKYP_10135.tif
## raster dimensions : 5860, 7806, 12 (nrow, ncol, nlyr)
## raster resolution : 500, 500 (x, y)
## raster extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## min values : 0.00645 0.00486 0.00368 0.00268 0.00203 0.00159 0.00128 0.00104 0.00084 0.00067 ...
## max values : 236 379 443 413 777 431 482 494 544 258 ...8.4.3 Plot all the area
# Plot Amazon in the year 2020
plot(soilMoisture.2020.rts, col=colorRamps::matlab.like2(1000))
# plot Amazon in October 2020.
plot(soilMoisture.2020.rts[['2020-10-01']], col=colorRamps::matlab.like2(1000), main="Soil Moisture in October 2020")
plot(amaz.basin.shp$geometry, add=T)
#
soilMoisture.2020.rts[['2020-10-01']] %>%
mask(mask = amaz.basin.shp) %>%
plot(col=colorRamps::matlab.like2(1000), main="Soil Moisture in October 2020")##
|---------|---------|---------|---------|
=========================================
plot(amaz.basin.shp$geometry, add=T)
8.4.4 Plot south of Amazon
x_min=-1.25e+06; x_max=0.05e+06; y_min=ymin(burntArea.rast); y_max=2.4e+06
# Plot south of Amazon in the year 2020.
soilMoisture.2020.rts@raster %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
mask(mask = amaz.basin.shp) %>%
plot(col=colorRamps::matlab.like2(1000))##
|---------|---------|---------|---------|
=========================================
|---------|---------|---------|---------|
=========================================
# Plot south of Amazon in October 2020.
soilMoisture.2020.rts[['2020-10-01']] %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
mask(mask = amaz.basin.shp) %>%
plot(col=colorRamps::matlab.like2(1000), main="Soil Moisture in October 2020")
plot(amaz.basin.shp$geometry, add=T)
8.5 Repalce negative values by NA for all data
# # Remove negative values.
# soilMoisture.rast[soilMoisture.rast < 0] <- NA
# # Save raster file
# writeRaster(soilMoisture.rast, "Amazon_new_data/4. Soil Moisture/soilmoisture_working_na.tif", overwrite=TRUE)
# Load data
soilMoisture.na.rast <- rast("Amazon_new_data/4. Soil Moisture/soilmoisture_working_na.tif")
soilMoisture.na.rast## class : SpatRaster
## dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## source : soilmoisture_working_na.tif
## names : soilm~01_01, soilm~01_02, soilm~01_03, soilm~01_04, soilm~01_05, soilm~01_06, ...
## min values : 3.900412e-02, 0.05741884, 0.07022829, 0.05813476, 0.1009746, 0.1354199, ...
## max values : 1.076058e+03, 750.97216797, 636.37585449, 721.50250244, 689.2701416, 540.9775391, ...8.6 Verification of the NA values
# cl <- makeCluster(detectCores() - 1)
# registerDoParallel(cl, cores=detectCores() - 1)
#
# tic("soilmoisture")
# #Raster to datatable in parallel: one raster per thread
# rasList <- foreach (ras_id=amaz.soilMoisture.list, .packages=c('terra', 'sf'), .combine='c') %dopar% {
# #Read all rasters into one big stack
# ras <- rast(ras_id)
# # Rename layers
# ras <- renameLayers(ras, 'soilmoisture_working_', '')
# # Replace negative value by `NA`
# ras[ras < 0] <- NA
# #
# ras.nonNA <- not.na(ras)
# ras.nonNA.mask <- mask(ras.nonNA, amaz.basin.shp)
# ras.freq.na <- freq(ras.nonNA.mask, digits=0, value=0, usenames=T)
#
# list(ras.freq.na)
# }
# stopCluster(cl)
# toc()
#
# #Bind all per-raster into one dataframe
# soilmoisture.freq.na <- rbindlist(rasList, fill=T, use.names=T)
# #
# colnames(soilmoisture.freq.na)[3] <- "soilmoisture_na"
# soilmoisture.freq.na <- soilmoisture.freq.na[order(soilmoisture.freq.na$layer)]
soilmoisture.freq.na8.7 Sample of NA values (February 2002)
8.7.1 Compute
soilmoisture.2002.02.name <- c("Amazon_new_data/4. Soil Moisture/03. Working Data/soilmoisture_working_2002_2.tif")
soilmoisture.2002.02.rast <- rast(soilmoisture.2002.02.name)
# Repalce negative values by `NA` for all data
soilmoisture.2002.02.rast[soilmoisture.2002.02.rast < 0] <- NA##
|---------|---------|---------|---------|
=========================================
|---------|---------|---------|---------|
=========================================
soilmoisture.2002.02.rast## class : SpatRaster
## dimensions : 5860, 7806, 1 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## source : spat_8KHrIMKbEZfP8DW_10135.tif
## name : soilm_2002_2_proj
## min value : 0.06081576
## max value : 465.46282959#
soilmoisture.2002.02.nonNA <- not.na(soilmoisture.2002.02.rast)##
|---------|---------|---------|---------|
=========================================
soilmoisture.2002.02.nonNA.mask <- mask(soilmoisture.2002.02.nonNA, amaz.basin.shp)##
|---------|---------|---------|---------|
=========================================
#
soilmoisture.2002.02.nonNA.df <- as.data.frame(soilmoisture.2002.02.nonNA.mask, xy=T)
soilmoisture.2002.02.NA.df <- soilmoisture.2002.02.nonNA.df[soilmoisture.2002.02.nonNA.df[,3]==F,]
soilmoisture.2002.02.NA.dfdim(soilmoisture.2002.02.NA.df)[1]## [1] 19758.7.2 Plot
# Plot
plot(soilmoisture.2002.02.nonNA.mask, col=c("darkorange1", "forestgreen"))
#
x_min=-0.45e+06; x_max=-0.3e+06; y_min=ymin(soilmoisture.2002.02.nonNA.mask); y_max=1.71e+06
soilmoisture.2002.02.nonNA.mask %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
plot(col=c("darkorange1", "forestgreen"))
#
x_min=1.2e+06; x_max=2.5e+06; y_min=3e+06; y_max=3.5e+06
soilmoisture.2002.02.nonNA.mask %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
plot(col=c("darkorange1", "forestgreen"))
#
x_min=0.1e+06; x_max=0.5e+06; y_min=4.25e+06; y_max=ymax(soilmoisture.2002.02.nonNA.mask)
soilmoisture.2002.02.nonNA.mask %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
plot(col=c("darkorange1", "forestgreen"))
#=====================
9 Elevation data
9.1 Import data
# list of files
amaz.elevation.list <- list.files("Amazon_new_data/5. Elevation/03. Working Data",
full.names=TRUE,
pattern = ".tif$")
# Import data with "Terra"
elevation.rast <- rast(amaz.elevation.list)
elevation.rast## class : SpatRaster
## dimensions : 5860, 7806, 1 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## source : elevation_working.tif
## name : elevation_proj
## min value : -85.5
## max value : 6470.59.2 Verification of the data
9.2.1 Verification the Names of layers
names(elevation.rast)## [1] "elevation_proj"9.2.2 Verification of the values
# Verification of the values
elevation.minmax <- minmax(elevation.rast) %>% t() %>% as.data.frame()
elevation.minmax9.3 Modification of the data
9.3.1 Rename layers
names(elevation.rast) <- "elevation"
elevation.rast## class : SpatRaster
## dimensions : 5860, 7806, 1 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## source : elevation_working.tif
## name : elevation
## min value : -85.5
## max value : 6470.59.3.2 Create ‘rts’ object
# Create a sequence date for 'rts' object
seq.dates <- seq(as.Date("2020-12-1"), as.Date("2020-12-1"), by = "month")
elevation.rts <- rts(elevation.rast, seq.dates)
elevation.rts## class : SpatRasterTS
## raster filename : /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/5. Elevation/03. Working Data/elevation_working.tif
## raster dimensions : 5860, 7806, 1 (nrow, ncol, nlyr)
## raster resolution : 500, 500 (x, y)
## raster extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## min values : -86
## max values : 64709.4 Plot
9.4.1 Plot all the area
# plot.
plot(elevation.rts[['2020-12-01']], col=colorRamps::matlab.like(1000), main="Elevation")
plot(amaz.basin.shp$geometry, add=T)
#
elevation.rts[['2020-12-01']] %>%
mask(mask = amaz.basin.shp) %>%
plot(col=colorRamps::matlab.like(1000), main="Elevation")##
|---------|---------|---------|---------|
=========================================
plot(amaz.basin.shp$geometry, add=T)
9.4.2 Plot south of Amazon
x_min=-1.25e+06; x_max=0.05e+06; y_min=ymin(burntArea.rast); y_max=2.4e+06
# Plot south of Amazon.
elevation.rts[['2020-12-01']] %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
mask(mask = amaz.basin.shp) %>%
plot(col=colorRamps::matlab.like(1000), main="Elevation")
plot(amaz.basin.shp$geometry, add=T)
9.5 Verification of the NA values
# cl <- makeCluster(detectCores() - 1)
# registerDoParallel(cl, cores=detectCores() - 1)
#
# tic("elevation")
# #Raster to datatable in parallel: one raster per thread
# rasList <- foreach (ras_id=amaz.elevation.list, .packages=c('terra', 'sf'), .combine='c') %dopar% {
# #Read all rasters into one big stack
# ras <- rast(ras_id)
# # Rename layers
# names(ras) <- "2020_12"
# #
# ras.nonNA <- not.na(ras)
# ras.nonNA.mask <- mask(ras.nonNA, amaz.basin.shp)
# ras.freq.na <- freq(ras.nonNA.mask, digits=0, value=0, usenames=T)
#
# list(ras.freq.na)
# }
# stopCluster(cl)
# toc()
#
# #Bind all per-raster into one dataframe
# elevation.freq.na <- rbindlist(rasList, fill=T, use.names=T)
# #
# colnames(elevation.freq.na)[3] <- "elevation_na"
# elevation.freq.na <- elevation.freq.na[order(elevation.freq.na$layer)]
elevation.freq.na#=====================
10 LandSurfaceTemp data
10.1 Import data
# list of files
amaz.landSurfaceTemp.list <- list.files("Amazon_new_data/6. LandSurfaceTemp/03. Working Data",
full.names=TRUE,
pattern = ".tif$")
# Import data with "Terra"
landSurfaceTemp.rast <- rast(amaz.landSurfaceTemp.list)
landSurfaceTemp.rast## class : SpatRaster
## dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## sources : landsurftemp_working_2001_1.tif
## landsurftemp_working_2001_10.tif
## landsurftemp_working_2001_11.tif
## ... and 237 more source(s)
## names : Month~ature, Month~ature, Month~ature, Month~ature, Month~ature, Month~ature, ...
## min values : 13255, 13746, 13376, 13395, 12869, 13508, ...
## max values : 16387, 16415, 16356, 16434, 16432, 16455, ...10.2 Verification of the data
10.2.1 Verification the Names of layers
names(landSurfaceTemp.rast) %>% head()## [1] "Monthly daytime 3min CMG Land-surface Temperature"
## [2] "Monthly daytime 3min CMG Land-surface Temperature"
## [3] "Monthly daytime 3min CMG Land-surface Temperature"
## [4] "Monthly daytime 3min CMG Land-surface Temperature"
## [5] "Monthly daytime 3min CMG Land-surface Temperature"
## [6] "Monthly daytime 3min CMG Land-surface Temperature"10.2.2 Verification of the values
# Verification of the values
landSurfaceTemp.minmax <- minmax(landSurfaceTemp.rast) %>% t() %>% as.data.frame()
landSurfaceTemp.minmax10.3 Modification of the data
10.3.1 Rename layers
landSurfaceTemp.rast <- renameLayers(landSurfaceTemp.rast, 'landsurftemp_working_', 'landsurftemp_')
# Order layers
ordered.names.6 <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month") %>%
format(., '%Y_%m') %>%
paste0("landsurftemp_", .)
landSurfaceTemp.rast <- landSurfaceTemp.rast[[ordered.names.6]]
# Print
landSurfaceTemp.rast## class : SpatRaster
## dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## sources : landsurftemp_working_2001_1.tif
## landsurftemp_working_2001_2.tif
## landsurftemp_working_2001_3.tif
## ... and 237 more source(s)
## names : lands~01_01, lands~01_02, lands~01_03, lands~01_04, lands~01_05, lands~01_06, ...
## min values : 13255, 12869, 13508, 13464, 13727, 13364, ...
## max values : 16387, 16432, 16455, 16223, 15990, 15897, ...cat('\n---------- \n\n')##
## ----------names(landSurfaceTemp.rast) %>% head()## [1] "landsurftemp_2001_01" "landsurftemp_2001_02" "landsurftemp_2001_03"
## [4] "landsurftemp_2001_04" "landsurftemp_2001_05" "landsurftemp_2001_06"10.3.2 Create ‘rts’ object
# Create a sequence date for 'rts' object
seq.dates <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month")
landSurfaceTemp.rts <- rts(landSurfaceTemp.rast, seq.dates)
landSurfaceTemp.rts## Raster Time Series with monthly periodicity from 2001-01-01 to 2020-12-01
## class : SpatRasterTS
## raster filename : /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/6. LandSurfaceTemp/03. Working Data/landsurftemp_working_2001_1.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/6. LandSurfaceTemp/03. Working Data/landsurftemp_working_2001_2.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/6. LandSurfaceTemp/03. Working Data/landsurftemp_working_2001_3.tif ...
## raster dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## raster resolution : 500, 500 (x, y)
## raster extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## min values : 13255 12869 13508 13464 13727 13364 13593 13570 13800 13746 ...
## max values : 16387 16432 16455 16223 15990 15897 15908 16023 16139 16415 ...10.4 Plot
10.4.1 Plot all the area
# Plot Amazon in the year 2020
plot(landSurfaceTemp.rts[['2020']], col=colorRamps::green2red(1000))
# plot Amazon in October 2020.
plot(landSurfaceTemp.rts[['2020-10-01']], col=colorRamps::green2red(1000), main="Land Surface Temp in October 2020")
plot(amaz.basin.shp$geometry, add=T)
#
landSurfaceTemp.rts[['2020-10-01']] %>%
mask(mask = amaz.basin.shp) %>%
plot(col=colorRamps::green2red(1000), main="Land Surface Temp in October 2020")##
|---------|---------|---------|---------|
=========================================
plot(amaz.basin.shp$geometry, add=T)
10.4.2 Plot south of Amazon
x_min=-1.25e+06; x_max=0.05e+06; y_min=ymin(burntArea.rast); y_max=2.4e+06
# Plot south of Amazon in the year 2020.
landSurfaceTemp.rts[['2020']]@raster %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
mask(mask = amaz.basin.shp) %>%
plot(col=colorRamps::green2red(1000))##
|---------|---------|---------|---------|
=========================================
|---------|---------|---------|---------|
=========================================
# Plot south of Amazon in October 2020.
landSurfaceTemp.rts[['2020-10-01']] %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
mask(mask = amaz.basin.shp) %>%
plot(col=colorRamps::green2red(1000), main="Land Surface Temp in October 2020")
plot(amaz.basin.shp$geometry, add=T)
10.5 Verification of the NA values
# cl <- makeCluster(detectCores() - 1)
# registerDoParallel(cl, cores=detectCores() - 1)
#
# tic("landsurftemp")
# #Raster to datatable in parallel: one raster per thread
# rasList <- foreach (ras_id=amaz.landSurfaceTemp.list, .packages=c('terra', 'sf'), .combine='c') %dopar% {
# #Read all rasters into one big stack
# ras <- rast(ras_id)
# # Rename layers
# ras <- renameLayers(ras, 'landsurftemp_working_', '')
# #
# ras.nonNA <- not.na(ras)
# ras.nonNA.mask <- mask(ras.nonNA, amaz.basin.shp)
# ras.freq.na <- freq(ras.nonNA.mask, digits=0, value=0, usenames=T)
#
# list(ras.freq.na)
# }
# stopCluster(cl)
# toc()
#
# #Bind all per-raster into one dataframe
# landsurftemp.freq.na <- rbindlist(rasList, fill=T, use.names=T)
# #
# colnames(landsurftemp.freq.na)[3] <- "landsurftemp_na"
# landsurftemp.freq.na <- landsurftemp.freq.na[order(landsurftemp.freq.na$layer)]
landsurftemp.freq.na10.6 Sample of NA values (February 2002)
10.6.1 Compute
landSurfaceTemp.2002.02.name <- c("Amazon_new_data/6. LandSurfaceTemp/03. Working Data/landsurftemp_working_2002_2.tif")
landSurfaceTemp.2002.02.rast <- rast(landSurfaceTemp.2002.02.name)
landSurfaceTemp.2002.02.rast## class : SpatRaster
## dimensions : 5860, 7806, 1 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## source : landsurftemp_working_2002_2.tif
## name : Monthly daytime 3min CMG Land-surface Temperature
## min value : 13421
## max value : 16484#
landSurfaceTemp.2002.02.nonNA <- not.na(landSurfaceTemp.2002.02.rast)##
|---------|---------|---------|---------|
=========================================
landSurfaceTemp.2002.02.nonNA.mask <- mask(landSurfaceTemp.2002.02.nonNA, amaz.basin.shp)##
|---------|---------|---------|---------|
=========================================
#
landSurfaceTemp.2002.02.nonNA.df <- as.data.frame(landSurfaceTemp.2002.02.nonNA.mask, xy=T)
landSurfaceTemp.2002.02.NA.df <- landSurfaceTemp.2002.02.nonNA.df[landSurfaceTemp.2002.02.nonNA.df[,3]==F,]
landSurfaceTemp.2002.02.NA.dfdim(landSurfaceTemp.2002.02.NA.df)[1]## [1] 208518810.6.2 Plot
# Plot
plot(landSurfaceTemp.2002.02.rast)
plot(landSurfaceTemp.2002.02.nonNA.mask, col=c("darkorange1", "forestgreen"))
plot(amaz.basin.shp$geometry, add=T)
#=====================
11 Specific Humidity data
11.1 Import data
# list of files
amaz.humidity.list <- list.files("Amazon_new_data/7. Specific Humidity/03. Working Data",
full.names=TRUE,
pattern = ".tif$")
# Import data with "Terra"
humidity.rast <- rast(amaz.humidity.list)
humidity.rast## class : SpatRaster
## dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## sources : humidity_working_2001_1.tif
## humidity_working_2001_10.tif
## humidity_working_2001_11.tif
## ... and 237 more source(s)
## names : humid~_proj, humid~_proj, humid~_proj, humid~_proj, humid~_proj, humid~_proj, ...
## min values : 0.003805317, 0.002230128, 0.002613974, 0.00270508, 0.004156957, 0.004153194, ...
## max values : 0.018998224, 0.019306520, 0.019978305, 0.02033914, 0.018878255, 0.019174447, ...11.2 Verification of the data
11.2.1 Verification the Names of layers
names(humidity.rast) %>% head()## [1] "humidity_2001_1_proj" "humidity_2001_10_proj" "humidity_2001_11_proj"
## [4] "humidity_2001_12_proj" "humidity_2001_2_proj" "humidity_2001_3_proj"11.2.2 Verification of the values
# Verification of the values
humidity.minmax <- minmax(humidity.rast) %>% t() %>% as.data.frame()
humidity.minmax11.3 Modification of the data
11.3.1 Rename layers
humidity.rast <- renameLayers(humidity.rast, 'humidity_working_', 'humidity_')
# Order layers
ordered.names.7 <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month") %>%
format(., '%Y_%m') %>%
paste0("humidity_", .)
humidity.rast <- humidity.rast[[ordered.names.7]]
# Print
humidity.rast## class : SpatRaster
## dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## sources : humidity_working_2001_1.tif
## humidity_working_2001_2.tif
## humidity_working_2001_3.tif
## ... and 237 more source(s)
## names : humid~01_01, humid~01_02, humid~01_03, humid~01_04, humid~01_05, humid~01_06, ...
## min values : 0.003805317, 0.004156957, 0.004153194, 0.003362617, 0.002170095, 0.001819314, ...
## max values : 0.018998224, 0.018878255, 0.019174447, 0.020020029, 0.019958658, 0.019542987, ...cat('\n---------- \n\n')##
## ----------names(humidity.rast) %>% head()## [1] "humidity_2001_01" "humidity_2001_02" "humidity_2001_03" "humidity_2001_04"
## [5] "humidity_2001_05" "humidity_2001_06"11.3.2 Create ‘rts’ object
# Create a sequence date for 'rts' object
seq.dates <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month")
humidity.rts <- rts(humidity.rast, seq.dates)
humidity.rts## Raster Time Series with monthly periodicity from 2001-01-01 to 2020-12-01
## class : SpatRasterTS
## raster filename : /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/7. Specific Humidity/03. Working Data/humidity_working_2001_1.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/7. Specific Humidity/03. Working Data/humidity_working_2001_2.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/7. Specific Humidity/03. Working Data/humidity_working_2001_3.tif ...
## raster dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## raster resolution : 500, 500 (x, y)
## raster extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## min values : 0.00381 0.00416 0.00415 0.00336 0.00217 0.00182 0.00111 0.00192 0.00218 0.00223 ...
## max values : 0.019 0.019 0.019 0.020 0.020 0.020 0.019 0.020 0.019 0.019 ...11.4 Plot
11.4.1 Plot all the area
# Plot Amazon in the year 2020
plot(humidity.rts[['2020']], col=topo.colors(20))
# plot Amazon in October 2020.
plot(humidity.rts[['2020-10-01']], col=topo.colors(20), main="Humidity in October 2020")
plot(amaz.basin.shp$geometry, add=T)
#
humidity.rts[['2020-10-01']] %>%
mask(mask = amaz.basin.shp) %>%
plot(col=topo.colors(20), main="Humidity in October 2020")
plot(amaz.basin.shp$geometry, add=T)
11.4.2 Plot south of Amazon
x_min=-1.25e+06; x_max=0.05e+06; y_min=ymin(burntArea.rast); y_max=2.4e+06
# Plot south of Amazon in the year 2020.
humidity.rts[['2020']]@raster %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
mask(mask = amaz.basin.shp) %>%
plot(col=topo.colors(20))##
|---------|---------|---------|---------|
=========================================
|---------|---------|---------|---------|
=========================================
# Plot south of Amazon in October 2020.
humidity.rts[['2020-10-01']] %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
mask(mask = amaz.basin.shp) %>%
plot(col=topo.colors(20), main="Humidity in October 2020")
plot(amaz.basin.shp$geometry, add=T)
11.5 Verification of the NA values
# cl <- makeCluster(detectCores() - 1)
# registerDoParallel(cl, cores=detectCores() - 1)
#
# tic("humidity")
# #Raster to datatable in parallel: one raster per thread
# rasList <- foreach (ras_id=amaz.humidity.list, .packages=c('terra', 'sf'), .combine='c') %dopar% {
# #Read all rasters into one big stack
# ras <- rast(ras_id)
# # Rename layers
# ras <- renameLayers(ras, 'humidity_working_', '')
# #
# ras.nonNA <- not.na(ras)
# ras.nonNA.mask <- mask(ras.nonNA, amaz.basin.shp)
# ras.freq.na <- freq(ras.nonNA.mask, digits=0, value=0, usenames=T)
#
# list(ras.freq.na)
# }
# stopCluster(cl)
# toc()
#
# #Bind all per-raster into one dataframe
# humidity.freq.na <- rbindlist(rasList, fill=T, use.names=T)
# #
# colnames(humidity.freq.na)[3] <- "humidity_na"
# humidity.freq.na <- humidity.freq.na[order(humidity.freq.na$layer)]
humidity.freq.na11.6 Sample of NA values (February 2002)
11.6.1 Compute
humidity.2002.02.name <- c("Amazon_new_data/7. Specific Humidity/03. Working Data/humidity_working_2002_2.tif")
humidity.2002.02.rast <- rast(humidity.2002.02.name)
humidity.2002.02.rast## class : SpatRaster
## dimensions : 5860, 7806, 1 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## source : humidity_working_2002_2.tif
## name : humidity_2002_2_proj
## min value : 0.004104336
## max value : 0.019729620#
humidity.2002.02.nonNA <- not.na(humidity.2002.02.rast)
humidity.2002.02.nonNA.mask <- mask(humidity.2002.02.nonNA, amaz.basin.shp)
#
humidity.2002.02.nonNA.df <- as.data.frame(humidity.2002.02.nonNA.mask, xy=T)
humidity.2002.02.NA.df <- humidity.2002.02.nonNA.df[humidity.2002.02.nonNA.df[,3]==F,]
humidity.2002.02.NA.dfdim(humidity.2002.02.NA.df)[1]## [1] 1816711.6.2 Plot
# Plot
plot(humidity.2002.02.nonNA.mask, col=c("darkorange1", "forestgreen"))
#
x_min=-0.45e+06; x_max=-0.3e+06; y_min=ymin(humidity.2002.02.nonNA.mask); y_max=1.71e+06
humidity.2002.02.nonNA.mask %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
plot(col=c("darkorange1", "forestgreen"))
#
x_min=0.e+06; x_max=xmax(humidity.2002.02.nonNA.mask); y_min=3.e+06; y_max=ymax(humidity.2002.02.nonNA.mask)
humidity.2002.02.nonNA.mask %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
plot(col=c("darkorange1", "forestgreen"))
#=====================
12 Evapotranspiration data
12.1 Import data
# list of files
amaz.evapotranspiration.list <- list.files("Amazon_new_data/8. Evapotranspiration/03. Working Data",
full.names=TRUE,
pattern = ".tif$")
# Import data with "Terra"
evapotranspiration.rast <- rast(amaz.evapotranspiration.list)
evapotranspiration.rast## class : SpatRaster
## dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## sources : evapotranspiration_working_2001_1.tif
## evapotranspiration_working_2001_10.tif
## evapotranspiration_working_2001_11.tif
## ... and 237 more source(s)
## names : evapo~_proj, evapo~_proj, evapo~_proj, evapo~_proj, evapo~_proj, evapo~_proj, ...
## min values : 0.000000e+00, -3.037881e-09, -1.953269e-08, -1.696159e-08, -3.509835e-08, -3.830653e-08, ...
## max values : 8.157189e-05, 8.477412e-05, 9.205872e-05, 7.918844e-05, 9.691325e-05, 8.340040e-05, ...12.2 Verification of the data
12.2.1 Verification the Names of layers
names(evapotranspiration.rast) %>% head()## [1] "evapotranspiration_2001_1_proj" "evapotranspiration_2001_10_proj"
## [3] "evapotranspiration_2001_11_proj" "evapotranspiration_2001_12_proj"
## [5] "evapotranspiration_2001_2_proj" "evapotranspiration_2001_3_proj"12.2.2 Verification of the values
# Verification of the values
evapotranspiration.minmax <- minmax(evapotranspiration.rast) %>% t() %>% as.data.frame()
evapotranspiration.minmax12.3 Modification of the data
12.3.1 Rename layers
evapotranspiration.rast <- renameLayers(evapotranspiration.rast, 'evapotranspiration_working_', 'evapotranspiration_')
# Order layers
ordered.names.8 <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month") %>%
format(., '%Y_%m') %>%
paste0("evapotranspiration_", .)
evapotranspiration.rast <- evapotranspiration.rast[[ordered.names.8]]
# Print
evapotranspiration.rast## class : SpatRaster
## dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## sources : evapotranspiration_working_2001_1.tif
## evapotranspiration_working_2001_2.tif
## evapotranspiration_working_2001_3.tif
## ... and 237 more source(s)
## names : evapo~01_01, evapo~01_02, evapo~01_03, evapo~01_04, evapo~01_05, evapo~01_06, ...
## min values : 0.000000e+00, -3.509835e-08, -3.830653e-08, -2.597946e-09, -1.043295e-08, -1.211072e-08, ...
## max values : 8.157189e-05, 9.691325e-05, 8.340040e-05, 7.563404e-05, 7.022559e-05, 8.013412e-05, ...cat('\n---------- \n\n')##
## ----------names(evapotranspiration.rast) %>% head()## [1] "evapotranspiration_2001_01" "evapotranspiration_2001_02"
## [3] "evapotranspiration_2001_03" "evapotranspiration_2001_04"
## [5] "evapotranspiration_2001_05" "evapotranspiration_2001_06"12.3.2 Create ‘rts’ object
# Create a sequence date for 'rts' object
seq.dates <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month")
evapotranspiration.rts <- rts(evapotranspiration.rast, seq.dates)
evapotranspiration.rts## Raster Time Series with monthly periodicity from 2001-01-01 to 2020-12-01
## class : SpatRasterTS
## raster filename : /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/8. Evapotranspiration/03. Working Data/evapotranspiration_working_2001_1.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/8. Evapotranspiration/03. Working Data/evapotranspiration_working_2001_2.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/8. Evapotranspiration/03. Working Data/evapotranspiration_working_2001_3.tif ...
## raster dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## raster resolution : 500, 500 (x, y)
## raster extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## min values : 0.0e+00 -3.5e-08 -3.8e-08 -2.6e-09 -1.0e-08 -1.2e-08 -5.4e-08 0.0e+00 -1.0e-08 -3.0e-09 ...
## max values : 8.2e-05 9.7e-05 8.3e-05 7.6e-05 7.0e-05 8.0e-05 8.3e-05 8.8e-05 7.9e-05 8.5e-05 ...12.4 Plot
12.4.1 Plot all the area
nbr_colors <- colorRampPalette(brewer.pal(9, "YlOrBr"))
# Plot Amazon in the year 2020
plot(evapotranspiration.rts[['2020']], col=nbr_colors(100))
# plot Amazon in October 2020.
plot(evapotranspiration.rts[['2020-10-01']], col=nbr_colors(100), main="Evapotranspiration in October 2020")
plot(amaz.basin.shp$geometry, add=T)
#
evapotranspiration.rts[['2020-10-01']] %>%
mask(mask = amaz.basin.shp) %>%
plot(col=nbr_colors(100), main="Evapotranspiration in October 2020")##
|---------|---------|---------|---------|
=========================================
plot(amaz.basin.shp$geometry, add=T)
12.4.2 Plot south of Amazon
nbr_colors <- colorRampPalette(brewer.pal(9, "YlOrBr"))
x_min=-1.25e+06; x_max=0.05e+06; y_min=ymin(burntArea.rast); y_max=2.4e+06
# Plot south of Amazon in the year 2020.
evapotranspiration.rts[['2020']]@raster %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
mask(mask = amaz.basin.shp) %>%
plot(col=nbr_colors(100))##
|---------|---------|---------|---------|
=========================================
|---------|---------|---------|---------|
=========================================
# Plot south of Amazon in October 2020.
evapotranspiration.rts[['2020-10-01']] %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
mask(mask = amaz.basin.shp) %>%
plot(col=nbr_colors(100), main="Evapotranspiration in October 2020")
plot(amaz.basin.shp$geometry, add=T)
12.5 Verification of the NA values
# cl <- makeCluster(detectCores() - 1)
# registerDoParallel(cl, cores=detectCores() - 1)
#
# tic("evapotranspiration")
# #Raster to datatable in parallel: one raster per thread
# rasList <- foreach (ras_id=amaz.evapotranspiration.list, .packages=c('terra', 'sf'), .combine='c') %dopar% {
# #Read all rasters into one big stack
# ras <- rast(ras_id)
# # Rename layers
# ras <- renameLayers(ras, 'evapotranspiration_working_', '')
# #
# ras.nonNA <- not.na(ras)
# ras.nonNA.mask <- mask(ras.nonNA, amaz.basin.shp)
# ras.freq.na <- freq(ras.nonNA.mask, digits=0, value=0, usenames=T)
#
# list(ras.freq.na)
# }
# stopCluster(cl)
# toc()
#
# #Bind all per-raster into one dataframe
# evapotranspiration.freq.na <- rbindlist(rasList, fill=T, use.names=T)
# #
# colnames(evapotranspiration.freq.na)[3] <- "evapotranspiration_na"
# evapotranspiration.freq.na <- evapotranspiration.freq.na[order(evapotranspiration.freq.na$layer)]
evapotranspiration.freq.na12.6 Sample of NA values (February 2002)
12.6.1 Compute
evapotranspi.2002.02.name <- c("Amazon_new_data/8. Evapotranspiration/03. Working Data/evapotranspiration_working_2002_2.tif")
evapotranspi.2002.02.rast <- rast(evapotranspi.2002.02.name)
evapotranspi.2002.02.rast## class : SpatRaster
## dimensions : 5860, 7806, 1 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## source : evapotranspiration_working_2002_2.tif
## name : evapotranspiration_2002_2_proj
## min value : -1.071572e-08
## max value : 7.709819e-05#
evapotranspi.2002.02.nonNA <- not.na(evapotranspi.2002.02.rast)##
|---------|---------|---------|---------|
=========================================
evapotranspi.2002.02.nonNA.mask <- mask(evapotranspi.2002.02.nonNA, amaz.basin.shp)##
|---------|---------|---------|---------|
=========================================
#
evapotranspi.2002.02.nonNA.df <- as.data.frame(evapotranspi.2002.02.nonNA.mask, xy=T)
evapotranspi.2002.02.NA.df <- evapotranspi.2002.02.nonNA.df[evapotranspi.2002.02.nonNA.df[,3]==F,]
evapotranspi.2002.02.NA.dfdim(evapotranspi.2002.02.NA.df)[1]## [1] 13299812.6.2 Plot
# Plot
plot(evapotranspi.2002.02.nonNA.mask, col=c("darkorange1", "forestgreen"))
#=====================
13 Wind Speed data
13.1 Import data
# list of files
amaz.wind.list <- list.files("Amazon_new_data/9. Wind Speed/03. Working Data",
full.names=TRUE,
pattern = ".tif$")
# Import data with "Terra"
wind.rast <- rast(amaz.wind.list)
wind.rast## class : SpatRaster
## dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## sources : wind_working_2001_1.tif
## wind_working_2001_10.tif
## wind_working_2001_11.tif
## ... and 237 more source(s)
## names : wind_~_proj, wind_~_proj, wind_~_proj, wind_~_proj, wind_~_proj, wind_~_proj, ...
## min values : 1.069407, 1.048391, 1.137132, 1.105208, 1.024364, 1.063414, ...
## max values : 8.202874, 7.747525, 7.947669, 8.901711, 9.390056, 8.901502, ...13.2 Verification of the data
13.2.1 Verification the Names of layers
names(wind.rast) %>% head()## [1] "wind_2001_1_proj" "wind_2001_10_proj" "wind_2001_11_proj"
## [4] "wind_2001_12_proj" "wind_2001_2_proj" "wind_2001_3_proj"13.2.2 Verification of the values
# Verification of the values
wind.minmax <- minmax(wind.rast) %>% t() %>% as.data.frame()
wind.minmax13.3 Modification of the data
13.3.1 Rename layers
wind.rast <- renameLayers(wind.rast, 'wind_working_', 'wind_')
# Order layers
ordered.names.9 <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month") %>%
format(., '%Y_%m') %>%
paste0("wind_", .)
wind.rast <- wind.rast[[ordered.names.9]]
# Print
wind.rast## class : SpatRaster
## dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## sources : wind_working_2001_1.tif
## wind_working_2001_2.tif
## wind_working_2001_3.tif
## ... and 237 more source(s)
## names : wind_2001_01, wind_2001_02, wind_2001_03, wind_2001_04, wind_2001_05, wind_2001_06, ...
## min values : 1.069407, 1.024364, 1.063414, 1.128978, 1.000217, 0.9819801, ...
## max values : 8.202874, 9.390056, 8.901502, 8.614891, 7.900918, 7.8636937, ...cat('\n---------- \n\n')##
## ----------names(wind.rast) %>% head()## [1] "wind_2001_01" "wind_2001_02" "wind_2001_03" "wind_2001_04" "wind_2001_05"
## [6] "wind_2001_06"13.3.2 Create ‘rts’ object
# Create a sequence date for 'rts' object
seq.dates <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month")
wind.rts <- rts(wind.rast, seq.dates)
wind.rts## Raster Time Series with monthly periodicity from 2001-01-01 to 2020-12-01
## class : SpatRasterTS
## raster filename : /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/9. Wind Speed/03. Working Data/wind_working_2001_1.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/9. Wind Speed/03. Working Data/wind_working_2001_2.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/9. Wind Speed/03. Working Data/wind_working_2001_3.tif ...
## raster dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## raster resolution : 500, 500 (x, y)
## raster extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## min values : 1.07 1.02 1.06 1.13 1.00 0.98 1.05 1.10 1.12 1.05 ...
## max values : 8.2 9.4 8.9 8.6 7.9 7.9 8.0 9.1 7.9 7.7 ...13.4 Plot
13.4.1 Plot all the area
# Plot Amazon in the year 2020
plot(wind.rts[['2020']], col=cm.colors(50))
# plot Amazon in October 2020.
plot(wind.rts[['2020-10-01']], col=cm.colors(50), main="Wind in October 2020")
plot(amaz.basin.shp$geometry, add=T)
#
wind.rts[['2020-10-01']] %>%
mask(mask = amaz.basin.shp) %>%
plot(col=cm.colors(50), main="Wind in October 2020")##
|---------|---------|---------|---------|
=========================================
plot(amaz.basin.shp$geometry, add=T)
13.4.2 Plot south of Amazon
x_min=-1.25e+06; x_max=0.05e+06; y_min=ymin(burntArea.rast); y_max=2.4e+06
# Plot south of Amazon in the year 2020.
wind.rts[['2020']]@raster %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
mask(mask = amaz.basin.shp) %>%
plot(col=cm.colors(50))##
|---------|---------|---------|---------|
=========================================
|---------|---------|---------|---------|
=========================================
# Plot south of Amazon in October 2020.
wind.rts[['2020-10-01']] %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
mask(mask = amaz.basin.shp) %>%
plot(col=cm.colors(50), main="Wind in October 2020")
plot(amaz.basin.shp$geometry, add=T)
13.5 Verification of the NA values
# cl <- makeCluster(detectCores() - 1)
# registerDoParallel(cl, cores=detectCores() - 1)
#
# tic("wind")
# #Raster to datatable in parallel: one raster per thread
# rasList <- foreach (ras_id=amaz.wind.list, .packages=c('terra', 'sf'), .combine='c') %dopar% {
# #Read all rasters into one big stack
# ras <- rast(ras_id)
# # Rename layers
# ras <- renameLayers(ras, 'wind_working_', '')
# #
# ras.nonNA <- not.na(ras)
# ras.nonNA.mask <- mask(ras.nonNA, amaz.basin.shp)
# ras.freq.na <- freq(ras.nonNA.mask, digits=0, value=0, usenames=T)
#
# list(ras.freq.na)
# }
# stopCluster(cl)
# toc()
#
# #Bind all per-raster into one dataframe
# wind.freq.na <- rbindlist(rasList, fill=T, use.names=T)
# #
# colnames(wind.freq.na)[3] <- "wind_na"
# wind.freq.na <- wind.freq.na[order(wind.freq.na$layer)]
wind.freq.na13.6 Sample of NA values (February 2002)
13.6.1 Compute
wind.2002.02.name <- c("Amazon_new_data/9. Wind Speed/03. Working Data/wind_working_2002_2.tif")
wind.2002.02.rast <- rast(wind.2002.02.name)
wind.2002.02.rast## class : SpatRaster
## dimensions : 5860, 7806, 1 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## source : wind_working_2002_2.tif
## name : wind_2002_2_proj
## min value : 1.199288
## max value : 9.328015#
wind.2002.02.nonNA <- not.na(wind.2002.02.rast)##
|---------|---------|---------|---------|
=========================================
wind.2002.02.nonNA.mask <- mask(wind.2002.02.nonNA, amaz.basin.shp)##
|---------|---------|---------|---------|
=========================================
#
wind.2002.02.nonNA.df <- as.data.frame(wind.2002.02.nonNA.mask, xy=T)
wind.2002.02.NA.df <- wind.2002.02.nonNA.df[wind.2002.02.nonNA.df[,3]==F,]
wind.2002.02.NA.dfdim(wind.2002.02.NA.df)[1]## [1] 1816713.6.2 Plot
# Plot
plot(wind.2002.02.nonNA.mask, col=c("darkorange1", "forestgreen"))
#
x_min=0.e+06; x_max=xmax(wind.2002.02.nonNA.mask); y_min=3.e+06; y_max=ymax(wind.2002.02.nonNA.mask)
wind.2002.02.nonNA.mask %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
plot(col=c("darkorange1", "forestgreen"))
#=====================
14 Air Temperature data
14.1 Import data
# list of files
amaz.airtemp.list <- list.files("Amazon_new_data/10. Air Temperature/03. Working Data",
full.names=TRUE,
pattern = ".tif$")
# Import data with "Terra"
airtemp.rast <- rast(amaz.airtemp.list)
airtemp.rast## class : SpatRaster
## dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## sources : airtemp_working_2001_1.tif
## airtemp_working_2001_10.tif
## airtemp_working_2001_11.tif
## ... and 237 more source(s)
## names : airte~_proj, airte~_proj, airte~_proj, airte~_proj, airte~_proj, airte~_proj, ...
## min values : 268.7742, 270.8836, 270.5267, 270.5608, 269.2101, 269.1326, ...
## max values : 303.1154, 304.1908, 304.2381, 303.4955, 304.7826, 304.7043, ...14.2 Verification of the data
14.2.1 Verification the Names of layers
names(airtemp.rast) %>% head()## [1] "airtemp_2001_1_proj" "airtemp_2001_10_proj" "airtemp_2001_11_proj"
## [4] "airtemp_2001_12_proj" "airtemp_2001_2_proj" "airtemp_2001_3_proj"14.2.2 Verification of the values
# Verification of the values
airtemp.minmax <- minmax(airtemp.rast) %>% t() %>% as.data.frame()
airtemp.minmax14.3 Modification of the data
14.3.1 Rename layers
airtemp.rast <- renameLayers(airtemp.rast, 'airtemp_working_', 'airtemp_')
# Order layers
ordered.names.10 <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month") %>%
format(., '%Y_%m') %>%
paste0("airtemp_", .)
airtemp.rast <- airtemp.rast[[ordered.names.10]]
# Print
airtemp.rast## class : SpatRaster
## dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## sources : airtemp_working_2001_1.tif
## airtemp_working_2001_2.tif
## airtemp_working_2001_3.tif
## ... and 237 more source(s)
## names : airte~01_01, airte~01_02, airte~01_03, airte~01_04, airte~01_05, airte~01_06, ...
## min values : 268.7742, 269.2101, 269.1326, 270.1584, 269.9435, 269.3073, ...
## max values : 303.1154, 304.7826, 304.7043, 304.4511, 303.5922, 304.0975, ...cat('\n---------- \n\n')##
## ----------names(airtemp.rast) %>% head()## [1] "airtemp_2001_01" "airtemp_2001_02" "airtemp_2001_03" "airtemp_2001_04"
## [5] "airtemp_2001_05" "airtemp_2001_06"14.3.2 Create ‘rts’ object
# Create a sequence date for 'rts' object
seq.dates <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month")
airtemp.rts <- rts(airtemp.rast, seq.dates)
airtemp.rts## Raster Time Series with monthly periodicity from 2001-01-01 to 2020-12-01
## class : SpatRasterTS
## raster filename : /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/10. Air Temperature/03. Working Data/airtemp_working_2001_1.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/10. Air Temperature/03. Working Data/airtemp_working_2001_2.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/10. Air Temperature/03. Working Data/airtemp_working_2001_3.tif ...
## raster dimensions : 5860, 7806, 240 (nrow, ncol, nlyr)
## raster resolution : 500, 500 (x, y)
## raster extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## min values : 269 269 269 270 270 269 269 269 270 271 ...
## max values : 303 305 305 304 304 304 304 305 305 304 ...14.4 Plot
14.4.1 Plot all the area
# Plot Amazon in the year 2020
plot(airtemp.rts[['2020']], col=colorRamps::matlab.like(100))
# plot Amazon in October 2020.
plot(airtemp.rts[['2020-10-01']], col=colorRamps::matlab.like(100), main="Air Temp. in October 2020")
plot(amaz.basin.shp$geometry, add=T)
#
airtemp.rts[['2020-10-01']] %>%
mask(mask = amaz.basin.shp) %>%
plot(col=colorRamps::matlab.like(100), main="Air Temp. in October 2020")##
|---------|---------|---------|---------|
=========================================
plot(amaz.basin.shp$geometry, add=T)
14.4.2 Plot south of Amazon
x_min=-1.25e+06; x_max=0.05e+06; y_min=ymin(burntArea.rast); y_max=2.4e+06
# Plot south of Amazon in the year 2020.
airtemp.rts[['2020']]@raster %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
mask(mask = amaz.basin.shp) %>%
plot(col=colorRamps::matlab.like(100))##
|---------|---------|---------|---------|
=========================================
|---------|---------|---------|---------|
=========================================
# Plot south of Amazon in October 2020.
airtemp.rts[['2020-10-01']] %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
mask(mask = amaz.basin.shp) %>%
plot(col=colorRamps::matlab.like(100), main="Air Temp. in October 2020")
plot(amaz.basin.shp$geometry, add=T)
14.5 Verification of the NA values
# cl <- makeCluster(detectCores() - 1)
# registerDoParallel(cl, cores=detectCores() - 1)
#
# tic("airtemp")
# #Raster to datatable in parallel: one raster per thread
# rasList <- foreach (ras_id=amaz.airtemp.list, .packages=c('terra', 'sf'), .combine='c') %dopar% {
# #Read all rasters into one big stack
# ras <- rast(ras_id)
# # Rename layers
# ras <- renameLayers(ras, 'airtemp_working_', '')
# #
# ras.nonNA <- not.na(ras)
# ras.nonNA.mask <- mask(ras.nonNA, amaz.basin.shp)
# ras.freq.na <- freq(ras.nonNA.mask, digits=0, value=0, usenames=T)
#
# list(ras.freq.na)
# }
# stopCluster(cl)
# toc()
#
# #Bind all per-raster into one dataframe
# airtemp.freq.na <- rbindlist(rasList, fill=T, use.names=T)
# #
# colnames(airtemp.freq.na)[3] <- "airtemp_na"
# airtemp.freq.na <- airtemp.freq.na[order(airtemp.freq.na$layer)]
airtemp.freq.na14.6 Sample of NA values (February 2002)
14.6.1 Compute
airtemp.2002.02.name <- c("Amazon_new_data/10. Air Temperature/03. Working Data/airtemp_working_2002_2.tif")
airtemp.2002.02.rast <- rast(airtemp.2002.02.name)
airtemp.2002.02.rast## class : SpatRaster
## dimensions : 5860, 7806, 1 (nrow, ncol, nlyr)
## resolution : 500, 500 (x, y)
## extent : -2156811, 1746189, 1625314, 4555314 (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic
## source : airtemp_working_2002_2.tif
## name : airtemp_2002_2_proj
## min value : 269.4591
## max value : 304.5671#
airtemp.2002.02.nonNA <- not.na(airtemp.2002.02.rast)##
|---------|---------|---------|---------|
=========================================
airtemp.2002.02.nonNA.mask <- mask(airtemp.2002.02.nonNA, amaz.basin.shp)##
|---------|---------|---------|---------|
=========================================
#
airtemp.2002.02.nonNA.df <- as.data.frame(airtemp.2002.02.nonNA.mask, xy=T)
airtemp.2002.02.NA.df <- airtemp.2002.02.nonNA.df[airtemp.2002.02.nonNA.df[,3]==F,]
airtemp.2002.02.NA.dfdim(airtemp.2002.02.NA.df)[1]## [1] 1816714.6.2 Plot
# Plot
plot(airtemp.2002.02.nonNA.mask, col=c("darkorange1", "forestgreen"))
#
x_min=0.e+06; x_max=xmax(airtemp.2002.02.nonNA.mask); y_min=3.e+06; y_max=ymax(airtemp.2002.02.nonNA.mask)
airtemp.2002.02.nonNA.mask %>%
crop(ext(x_min, x_max, y_min, y_max)) %>%
plot(col=c("darkorange1", "forestgreen"))
#=====================
15 Merge all dataframes of NA values.
# create the dataframe
seq.dates.str <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month") %>% format("%Y_%m")
amaz.na.df <- as.data.frame(seq.dates.str)
colnames(amaz.na.df) <- "layer"
# Merge
amaz.na.df <- list(amaz.na.df,
burntArea.freq.na[,-2],
landCover.freq.na[,-2],
precipitation.freq.na[,-2],
soilmoisture.freq.na[,-2],
elevation.freq.na[,-2],
landsurftemp.freq.na[,-2],
humidity.freq.na[,-2],
evapotranspiration.freq.na[,-2],
wind.freq.na[,-2],
airtemp.freq.na[,-2]) %>%
reduce(full_join, by="layer")
amaz.na.dfamaz.na.df[133:144,]#=====================